{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 第15章 奇异值分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 习题15.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "&emsp;&emsp;试求矩阵\n",
    "$$\n",
    "A = \\left[ \n",
    "\\begin{array}{lll}\n",
    "1 & 2 & 0 \\\\\n",
    "2 & 0 & 2 \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$的奇异值分解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**\n",
    "\n",
    "1. 给出奇异值分解的定义\n",
    "2. 给出奇异值分解的步骤\n",
    "3. 使用`numpy`实现奇异值分解\n",
    "4. 自编程实现奇异值分解\n",
    "5. 推导计算奇异值分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：奇异值分解**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第15.1节的奇异值分解的定义：\n",
    "\n",
    "> **定义15.1（奇异值分解）** 矩阵的奇异值分解是指，将一个非零的$m \\times n$实矩阵$A$，$A \\in R^{m \\times n}$，表示为以下三个实矩阵乘积形式的运算，即进行矩阵的因子分解：\n",
    "> $$ A = U \\Sigma V^T $$\n",
    "> 其中$U$是$m$阶正交矩阵，$V$是$n$阶正交矩阵，$\\Sigma$是由降序排列的非负的对角线元素组成的$m \\times n$矩形对角矩阵，\n",
    "> $$ \\begin{array}{l}\n",
    "U U^T = I \\\\\n",
    "V V^T = I \\\\\n",
    "\\Sigma = \\text{diag}(\\sigma_1, \\sigma_2, \\cdots, \\sigma_p) \\\\\n",
    "\\sigma_1 \\geqslant \\sigma_2 \\geqslant \\cdots \\sigma_p \\geqslant 0 \\\\\n",
    "p = \\min(m, n)\n",
    "\\end{array} $$\n",
    "> $U \\Sigma V^T$称为矩阵$A$的奇异值分解，$\\sigma_i$称为矩阵$A$的奇异值，$U$的列向量称为左奇异向量，$V$的列向量称为右奇异向量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：奇异值分解的一般步骤**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第15.2节的奇异值分解的计算：\n",
    "\n",
    "> 给定$m \\times n$矩阵$A$，可以写出矩阵奇异值分解的计算过程。  \n",
    "> （1）求$A^T A$的特征值和特征向量。  \n",
    "> 计算对称矩阵$W=A^T A$。  \n",
    "> 求解特征方程\n",
    "> $$ (W - \\lambda I)x = 0 $$\n",
    "> 得到特征值$\\lambda_i$，并将特征值由大到小排列\n",
    "> $$ \\lambda_1 \\geqslant \\lambda_2 \\geqslant \\cdots \\geqslant \\lambda_n \\geqslant 0 $$\n",
    "> 将特征值$\\lambda_i \\ (i=1,2,\\cdots,n)$代入特征方程求得对应的特征向量。  \n",
    "> （2）求$n$阶正交矩阵$V$  \n",
    "> 将特征向量单位化，得到单位特征向量$v_1, v_2, \\cdots, v_n$，构成$n$阶正交矩阵$V$  \n",
    "> $$ V = [v_1 \\ v_2 \\ \\cdots \\ v_n] $$\n",
    "> （3）求$m \\times n$对角矩阵$\\Sigma$  \n",
    "> 计算$A$的奇异值\n",
    "> $$ \\sigma_i = \\sqrt{\\lambda_i}, \\quad i=1,2,\\cdots,n $$\n",
    "> 构造$m \\times n$矩形对角矩阵$\\Sigma$，主对角线元素是奇异值，其余元素是零\n",
    "> $$ \\Sigma = \\text{diag} (\\sigma_1, \\sigma_2, \\cdots, \\sigma_n) $$\n",
    "> （4）求$m$阶正交矩阵$U$  \n",
    "> 对$A$的前$r$个正奇异值，令  \n",
    "> $$ u_j = \\frac{1}{\\sigma_j} A v_j, \\quad j = 1, 2, \\cdots, r $$\n",
    "> 得到\n",
    "> $$ U_1 = [u1 \\ u2 \\ \\cdots \\ u_r] $$\n",
    "> 求$A^T$的零空间的一组标准正交基$\\{u_{r+1}, u_{r+2}, \\cdots, u_m \\}$，令\n",
    "> $$ U_2 = [u_{r+1} \\ u_{r+2} \\ \\cdots \\ u_m ] $$\n",
    "> 并令\n",
    "> $$ U = [U_1 \\ U_2] $$\n",
    "> （5）得到奇异值分解\n",
    "> $$ A = U \\Sigma V^T $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：使用numpy实现奇异值分解**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[ 0.45  0.89]\n",
      " [ 0.89 -0.45]]\n",
      "S= [3. 2.]\n",
      "V= [[ 0.75 -0.   -0.67]\n",
      " [ 0.3   0.89  0.33]\n",
      " [ 0.6  -0.45  0.67]]\n",
      "A= [[ 1.  2.  0.]\n",
      " [ 2. -0.  2.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.array([[1, 2, 0],\n",
    "              [2, 0, 2]])\n",
    "\n",
    "# 调用numpy的svd方法\n",
    "U, S, V = np.linalg.svd(A)\n",
    "\n",
    "# 设置精度为2\n",
    "np.set_printoptions(precision=2, suppress=True)\n",
    "\n",
    "print(\"U=\", U)\n",
    "print(\"S=\", S)\n",
    "print(\"V=\", V.T)\n",
    "Sigma = np.zeros_like(A, float)\n",
    "np.fill_diagonal(Sigma, S)\n",
    "calc = np.dot(np.dot(U, Sigma), V)\n",
    "print(\"A=\", calc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第4步：自编程实现奇异值分解**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import null_space\n",
    "\n",
    "\n",
    "def my_svd(A):\n",
    "    m = A.shape[0]\n",
    "\n",
    "    # (1) 计算对称矩阵 A^T A 的特征值与特征向量，\n",
    "    W = np.dot(A.T, A)\n",
    "    # 返回的特征值lambda_value是升序的，特征向量V是单位化的特征向量\n",
    "    lambda_value, V = np.linalg.eigh(W)\n",
    "    # 并按特征值从大到小排列\n",
    "    lambda_value = lambda_value[::-1]\n",
    "    lambda_value = lambda_value[lambda_value > 0]\n",
    "    # (2)计算n阶正价矩阵V\n",
    "    V = V[:, -1::-1]\n",
    "\n",
    "    # (3) 求 m * n 对角矩阵\n",
    "    sigma = np.sqrt(lambda_value)\n",
    "    S = np.diag(sigma) @ np.eye(*A.shape)\n",
    "\n",
    "    # (4.1) 求A的前r个正奇异值\n",
    "    r = np.linalg.matrix_rank(A)\n",
    "    U1 = np.hstack([(np.dot(A, V[:, i]) / sigma[i])[:, np.newaxis] for i in range(r)])\n",
    "    # (4.2) 求A^T的零空间的一组标准正交基\n",
    "    U = U1\n",
    "    if r < m:\n",
    "        U2 = null_space(A.T)\n",
    "        U2 = U2[:, r:]\n",
    "        U = np.hstack([U, U2])\n",
    "\n",
    "    return U, S, V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[-0.45 -0.89]\n",
      " [-0.89  0.45]]\n",
      "S= [[3. 0. 0.]\n",
      " [0. 2. 0.]]\n",
      "V= [[-0.75 -0.   -0.67]\n",
      " [-0.3  -0.89  0.33]\n",
      " [-0.6   0.45  0.67]]\n",
      "A= [[ 1.  2. -0.]\n",
      " [ 2. -0.  2.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 2, 0],\n",
    "              [2, 0, 2]])\n",
    "\n",
    "np.set_printoptions(precision=2, suppress=True)\n",
    "\n",
    "U, S, V = my_svd(A)\n",
    "print(\"U=\", U)\n",
    "print(\"S=\", S)\n",
    "print(\"V=\", V)\n",
    "calc = np.dot(np.dot(U, S), V.T)\n",
    "print(\"A=\", calc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第5步：推导计算奇异值分解**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 求解$A^TA$的特征值和特征向量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;计算对称矩阵$W = A^TA$\n",
    "\n",
    "$$\n",
    "A^TA = \n",
    "\\left[\\begin{array}{cc}\n",
    "1 & 2 \\\\\n",
    "2 & 0 \\\\\n",
    "0 & 2\n",
    "\\end{array}\\right]\n",
    "\\left[\\begin{array}{ccc}\n",
    "1 & 2 & 0 \\\\\n",
    "2 & 0 &2 \n",
    "\\end{array}\\right]\n",
    "=\n",
    "\\left[\\begin{array}{ccc}\n",
    "5 & 2 & 4 \\\\\n",
    "2 & 4 & 0 \\\\\n",
    "4 & 0 & 4\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;特征值和特征向量满足方程\n",
    "\n",
    "$$\n",
    "( A^TA- \\lambda I)x=0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;得到齐次方程组\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{ccc}\n",
    "5-\\lambda & 2 & 4 \\\\\n",
    "2 & 4-\\lambda & 0 \\\\\n",
    "4 & 0 & 4-\\lambda\n",
    "\\end{array}\\right]x = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;该方程组有非零解的充要条件是\n",
    "\n",
    "$$\n",
    "\\left| \n",
    "\\begin{array}{ccc}\n",
    "5-\\lambda & 2 & 4 \\\\\n",
    "2 & 4-\\lambda & 0 \\\\\n",
    "4 & 0 & 4-\\lambda\n",
    "\\end{array} \\right| = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;解得矩阵 $W = A^TA$ 的特征值 $\\lambda_1=9$、$\\lambda_2=4$ 和 $\\lambda_3=0$ （按从大到小排列）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;特征值代入对应的齐次方程组，得到对应的单位特征向量 **(特征向量的选取不唯一）**："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "v_1 = \\left[ \\begin{array}{c}\n",
    "\\displaystyle \\frac{5}{3\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{2}{3\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{4}{3\\sqrt{5}}\n",
    "\\end{array} \\right],\n",
    "v_2 = \\left[ \\begin{array}{c}\n",
    "0 \\\\\n",
    "\\displaystyle -\\frac{2}{\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{1}{\\sqrt{5}}\n",
    "\\end{array} \\right],\n",
    "v_3 = \\left[ \\begin{array}{c}\n",
    "\\displaystyle -\\frac{2}{3}\\\\\n",
    "\\displaystyle \\frac{1}{3} \\\\\n",
    "\\displaystyle \\frac{2}{3}\n",
    "\\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. 求正交矩阵$V$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;构造正交矩阵$V$\n",
    "$$\n",
    "V = \\left[ \\begin{array}{ccc}\n",
    "v_1 & v_2 & v_3 \n",
    "\\end{array} \\right] \n",
    "=\n",
    "\\left[ \\begin{array}{ccc}\n",
    "\\displaystyle \\frac{5}{3\\sqrt{5}} & 0 & \\displaystyle -\\frac{2}{3} \\\\\n",
    "\\displaystyle \\frac{2}{3\\sqrt{5}} & \\displaystyle -\\frac{2}{\\sqrt{5}} & \\displaystyle \\frac{1}{3} \\\\\n",
    "\\displaystyle \\frac{4}{3\\sqrt{5}} & \\displaystyle \\frac{1}{\\sqrt{5}} & \\displaystyle \\frac{2}{3}\n",
    "\\displaystyle \\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. 求解对角矩阵$\\Sigma$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;构造对角矩阵$\\Sigma$\n",
    "$$\n",
    "\\Sigma = \n",
    "\\left[ \\begin{array}{ccc}\n",
    "3 & 0 & 0 \\\\\n",
    "0 & 2 & 0 \n",
    "\\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. 求正交矩阵$U$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;基于$A$的正奇异值 $\\sigma_1=3$ 、$\\sigma_2=2$ 计算得到列向量$u_1,u_2$：\n",
    "\n",
    "$$\n",
    "u_1= \\frac{1}{\\sigma_1}A v_1 \n",
    "= \\frac{1}{3}\n",
    "\\left[ \\begin{array}{ccc}\n",
    "1 & 2 & 0 \\\\\n",
    "2 & 0 &2 \n",
    "\\end{array} \\right]\n",
    "\\left[ \\begin{array}{c}\n",
    "\\displaystyle \\frac{5}{3\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{2}{3\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{4}{3\\sqrt{5}}\n",
    "\\end{array} \\right] \n",
    "= \\left[ \\begin{array}{c}\n",
    "\\displaystyle \\frac{1}{\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{2}{\\sqrt{5}} \n",
    "\\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u_2 = \\frac{1}{\\sigma_2} A v_2 = \n",
    "\\frac{1}{2}\n",
    "\\left[ \\begin{array}{ccc}\n",
    "1 & 2 & 0 \\\\\n",
    "2 & 0 &2 \n",
    "\\end{array} \\right]\n",
    "\\left[ \\begin{array}{c}\n",
    "0 \\\\\n",
    "\\displaystyle -\\frac{2}{\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{1}{\\sqrt{5}}\n",
    "\\end{array} \\right] \n",
    "= \\left[ \\begin{array}{c}\n",
    "\\displaystyle -\\frac{2}{\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{1}{\\sqrt{5}} \n",
    "\\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;构造正交矩阵$U$\n",
    "\n",
    "$$\n",
    "U= \\left[ \\begin{array}{cc}\n",
    "\\displaystyle \\frac{1}{\\sqrt{5}} & \\displaystyle -\\frac{2}{\\sqrt{5}}\\\\\n",
    "\\displaystyle \\frac{2}{\\sqrt{5}} & \\displaystyle \\frac{1}{\\sqrt{5}} \n",
    "\\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. 矩阵$A$的奇异值分解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A = U \\Sigma V^T = \n",
    "\\left[ \\begin{array}{cc}\n",
    "\\displaystyle \\frac{1}{\\sqrt{5}} & \\displaystyle -\\frac{2}{\\sqrt{5}} \\\\\n",
    "\\displaystyle \\frac{2}{\\sqrt{5}} & \\displaystyle \\frac{1}{\\sqrt{5}} \n",
    "\\end{array} \\right]\n",
    "\\left[ \\begin{array}{ccc}\n",
    "3 & 0 & 0 \\\\\n",
    "0 & 2 & 0 \n",
    "\\end{array} \\right]\n",
    "\\left[ \\begin{array}{ccc}\n",
    "\\displaystyle \\frac{5}{3\\sqrt{5}} & 0 & \\displaystyle -\\frac{2}{3} \\\\\n",
    "\\displaystyle \\frac{2}{3\\sqrt{5}} & \\displaystyle -\\frac{2}{\\sqrt{5}} & \\displaystyle \\frac{1}{3} \\\\\n",
    "\\displaystyle \\frac{4}{3\\sqrt{5}} & \\displaystyle \\frac{1}{\\sqrt{5}} & \\displaystyle \\frac{2}{3}\n",
    "\\end{array} \\right]^T\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题15.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;试求矩阵\n",
    "$$\n",
    "A = \\left[ \n",
    "\\begin{array}{lll}\n",
    "2 & 4 \\\\\n",
    "1 & 3 \\\\\n",
    "0 & 0 \\\\\n",
    "0 & 0 \n",
    "\\end{array}\n",
    "\\right]\n",
    "$$的奇异值分解并写出其外积展开式。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**  \n",
    "\n",
    "1. 给出矩阵的外积展开式定义\n",
    "2. 使用`numpy`计算矩阵$A$的奇异值分解\n",
    "3. 根据奇异值分解的结果，写出外积展开式 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：矩阵 $A$ 的外积展开式；**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第15.3.3节的矩阵的外积展开式：\n",
    "\n",
    "> &emsp;&emsp;矩阵$A$的奇异值分解$U \\Sigma V^T$也可以由外积形式表示。事实上，若将$A$的奇异值分解看成矩阵$U \\Sigma$和$V^T$的乘积，将$U \\Sigma$按列向量分块，将$V^T$按行向量分块，即得\n",
    "> $$ \\begin{array}{l}\n",
    "U \\Sigma = \\big[ \\sigma_1 u_1 \\ \\sigma_2 u_2 \\ \\cdots \\ \\sigma_n u_n \\big] \\\\\n",
    "V^T = \\left[ \\begin{array}{c} \n",
    "v_1^T \\\\\n",
    "v_2^T \\\\\n",
    "\\vdots \\\\\n",
    "v_n^T \n",
    "\\end{array} \\right]\n",
    "\\end{array} $$\n",
    "> 则\n",
    "> $$ A = \\sigma_1 u_1 v_1^T + \\sigma_2 u_2 v_2^T + \\cdots + \\sigma_n u_n v_n^T $$\n",
    "> 称为矩阵$A$的外积展开式，其中$u_k v_k^T$为$m \\times n$矩阵，是列向量$u_k$和行向量$v_k^T$的外积，其第$i$行第$j$列元素为$u_k$的第$i$个元素与$v_k^T$的第$j$个元素的乘积。即\n",
    "> $$ u_i v_j^T = \\left[ \n",
    "\\begin{array}{c} \n",
    "u_{1i} \\\\\n",
    "u_{2i} \\\\\n",
    "\\vdots \\\\\n",
    "u_{mi} \n",
    "\\end{array} \n",
    "\\right] \\big[ v_{1j} \\ v_{2j} \\ \\cdots \\ v_{nj} \\big] \n",
    "= \\left[ \\begin{array}{cccc} \n",
    "u_{1i} v_{1j} & u_{1i} v_{2j} & \\cdots & u_{1i} v_{nj} \\\\\n",
    "u_{2i} v_{1j} & u_{2i} v_{2j} & \\cdots & u_{2i} v_{nj} \\\\\n",
    "\\vdots & \\vdots & & \\vdots \\\\\n",
    "u_{mi} v_{1j} & u_{mi} v_{2j} & \\cdots & u_{mi} v_{nj}\n",
    "\\end{array} \\right] $$\n",
    "> $A$的外积展开式也可以写成下面的形式\n",
    "> $$ A = \\sum_{k=1}^n A_k = \\sum_{k=1}^n \\sigma_k u_k v_k^T $$\n",
    "> 其中$A_k = \\sigma_k u_k v_k^T$是$m \\times n$矩阵。上式将矩阵$A$分解为矩阵的有序加权和。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第14章的本章概要中的外积展开式：\n",
    "\n",
    "> 任意一个实矩阵$A$可以由其外积展开式表示\n",
    "> $$ A = \\sigma_1 u_1 v_1^T + \\sigma_2 u_2 v_2^T + \\cdots + \\sigma_n u_n v_n^T $$\n",
    "> 其中$u_k v_k^T$为$m \\times n$矩阵，是列向量$u_k$和行向量$v_k^T$的外积，$\\sigma_k$为奇异值，$u_k, v_k^T, \\sigma_k$通过矩阵$A$的奇异值分解得到。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：计算矩阵$A$的奇异值分解**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[-0.82 -0.58  0.    0.  ]\n",
      " [-0.58  0.82  0.    0.  ]\n",
      " [ 0.    0.    1.    0.  ]\n",
      " [ 0.    0.    0.    1.  ]]\n",
      "S= [5.46 0.37]\n",
      "V= [[-0.4  -0.91]\n",
      " [-0.91  0.4 ]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.array([[2, 4],\n",
    "              [1, 3],\n",
    "              [0, 0],\n",
    "              [0, 0]])\n",
    "\n",
    "# 调用numpy的svd方法\n",
    "U, S, V = np.linalg.svd(A)\n",
    "np.set_printoptions()\n",
    "\n",
    "print(\"U=\", U)\n",
    "print(\"S=\", S)\n",
    "print(\"V=\", V.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：根据奇异值分解的结果，写出外积展开式**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;矩阵$A$的外积展开式为\n",
    "$$\n",
    "A = \\sigma_1 u_1 v_1^T + \\sigma_2 u_2 v_2^T\n",
    "$$ \n",
    "其中\n",
    "$$\n",
    "\\sigma_1 = 5.4649857, \\sigma_2 = 0.36596619 \\\\ \n",
    "u_1 = \\left [ \\begin{array}{c}\n",
    "-0.81741556 \\\\ \n",
    "-0.57604844 \\\\\n",
    "0 \\\\\n",
    "0\n",
    "\\end{array} \\right],\n",
    "u_2 = \\left[ \\begin{array}{c}\n",
    "-0.57604844 \\\\\n",
    "0.81741556 \\\\\n",
    "0 \\\\\n",
    "0\n",
    "\\end{array} \\right ] \\\\\n",
    "v_1^T = [ -0.40455358, -0.9145143 ] , v_2^T = [ -0.9145143,0.40455358 ]\n",
    "$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A= [[ 2.  4.]\n",
      " [ 1.  3.]\n",
      " [-0.  0.]\n",
      " [-0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "calc = S[0] * np.outer(U[:, 0], V[:, 0]) + S[1] * np.outer(U[:, 1], V[:, 1])\n",
    "print(\"A=\", calc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题15.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;比较矩阵的奇异值分解与对称矩阵的对角化的异同。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**  \n",
    "\n",
    "1. 给出矩阵的奇异值分解的定义\n",
    "2. 给出对称矩阵和方阵的对角化定义\n",
    "3. 比较两者的相同点\n",
    "4. 比较两者的不同点"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**注意：** 与书中一致，不考虑复数矩阵。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：矩阵的奇异值分解的定义：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;矩阵的奇异值分解定义，详见书中第271页的定义15.1（奇异值分解）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：对称矩阵和方阵对角化定义；**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据《实对称矩阵对角化教学的应用案例》（张丽静、刘白羽、申亚男，2019）对称矩阵的对角化：\n",
    "\n",
    "> 设$A$为$n$阶实对称矩阵，则存在一个正交矩阵$P = (p_1, p_2, \\cdots, p_n)$，使得\n",
    "> \n",
    "> $$ A = P \\Lambda P^T $$\n",
    "> \n",
    "> 其中，$\\Lambda$为对角元素是$A$的特征值$\\lambda_1,\\lambda_2,\\cdots,\\lambda_n$所构成的对角矩阵；$p_1,p_2,\\cdots,p_n$为特征值$\\lambda_1,\\lambda_2,\\cdots,\\lambda_n$对应的单位特征向量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据参考资料：https://www.matongxue.com/parts/4628 中关于方阵的对角化：\n",
    "> 如果$n$阶方阵$A$有$n$个线性无关的特征向量$p_1, p_2, \\cdots,p_n$，那么如下矩阵：\n",
    "> \n",
    "> $$ P= (p_1, p_2, \\cdots, p_n) $$\n",
    "> \n",
    "> 可以使得：$ A = P \\Lambda P^{-1}$，其中$\\Lambda$为对角阵$\\Lambda = \\text{diag} \\{ \\lambda_1, \\lambda_2, \\cdots, \\lambda_n \\}$，特征值 $\\lambda_1,\\lambda_2,\\cdots, \\lambda_n$为特征向量$p_1, p_2, \\cdots, p_n$对应的特征值，该过程为对角化（Diagonalizable）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：两者的相同点**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 两者都可以实现特征分解，都需要求解特征值与特征向量\n",
    "2. 实对称矩阵的对角化是方阵对角化的特例，矩阵奇异值分解又是方阵对角化的推广\n",
    "3. 实对称矩阵一定可以对角化，矩阵的奇异值分解也一定存在\n",
    "4. 都可以被正交矩阵对角化\n",
    "5. 本质上都是用不同的基来描述线性变化"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第4步：两者的不同点**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 矩阵对角化要求矩阵$A$是一个方阵，矩阵奇异值分解不需要对矩阵$A$有特殊要求\n",
    "2. 奇异值分解中要求对角元素降序排列，对称矩阵的对角化没有这一要求"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题15.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;证明任何一个秩为1的矩阵可以写成两个向量的外积形式，并给出实例。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 假设矩阵$A$的秩为1，写出其奇异值分解形式\n",
    "2. 将奇异值分解的矩阵$\\Sigma$写成外积形式\n",
    "3. 将矩阵$A$的奇异值分解写成外积形式\n",
    "4. 给出实例"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：假设矩阵 $A$ 的秩为1，写出其奇异值分解形式**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;设矩阵$A_{m\\times n}$的秩为1。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据书中第15章的定理15.1的奇异值分解基本定理：\n",
    "\n",
    "> &emsp;&emsp;**定理15.1（奇异值分解基本定理）** 若$A$为一$m \\times n$实矩阵，$A \\in R^{m \\times n}$，则$A$的奇异值分解存在\n",
    "> $$ A = U \\Sigma V^T $$\n",
    "> 其中$U$是$m$阶正交矩阵，$V$是$n$阶正交矩阵，$\\Sigma$是$m \\times n$矩形对角矩阵，其对角线元素非负，且按降序排列。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;则存在矩阵$A$的奇异值分解\n",
    "$$\n",
    "A= U \\Sigma V^T \\tag{1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：将奇异值分解的矩阵$\\Sigma$写成外积形式；** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据矩阵$A$的秩与对角矩阵$\\Sigma$的秩之间的关系：\n",
    "$$\n",
    "\\text{rank}(\\Sigma) = \\text{rank}(A) = 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据奇异值分解基本定理，对角矩阵$\\Sigma$的奇异值是降序排列的，所以其中的非零元素一定位于第一行第一列，可设为 \n",
    "\n",
    "$$\n",
    "\\left[ \\begin{array}{cc}\n",
    "\\sigma_1 & \\\\ \n",
    "& O_{(m-1) \\times (n-1)} \n",
    "\\end{array} \\right]_{m \\times n}\n",
    "$$\n",
    "\n",
    "其中，$O_{(m-1) \\times (n-1)}$表示$m-1$行$n-1$列的零矩阵。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;设两向量 $a = [1,0,\\cdots,0]_{m\\times 1}^T, b = [\\sigma_1, 0, \\cdots, 0]_{n\\times 1}^T$，可得：\n",
    "\n",
    "$$\n",
    "\\Sigma = a b^T \\tag{2}\n",
    "$$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：将矩阵$A$ 的奇异值分解写成外积形式** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;将式（2）代入式（1）中可得：\n",
    "\n",
    "$$\n",
    "A = U a b^T V^T = (Ua)(Vb)^T \n",
    "$$ \n",
    "\n",
    "其中，$Ua$是$m\\times 1$阶的列向量，$Vb$是$n\\times 1$阶的列向量。\n",
    "\n",
    "&emsp;&emsp;所以，矩阵$A$可以写成向量$Ua$和向量$Vb$的外积形式，命题得证。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第4步：给出实例** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;使用书中第283页例15.5的矩阵$A = \\left[ \\begin{array}{ccc} 1 & 1 \\\\ 2 & 2 \\\\ 0 & 0 \\end{array} \\right]$的秩为1，矩阵$A$的奇异值分解为\n",
    "$$\n",
    "\\begin{aligned}\n",
    "A \n",
    "&= U\\Sigma V^T \\\\\n",
    "&= \\left[ \\begin{array}{ccc}\\frac{1}{\\sqrt{5}}&-\\frac{2}{\\sqrt{5}}&0\\\\\\frac{2}{\\sqrt{5}}&\\frac{1}{\\sqrt{5}}&0\\\\0&0&1\\end{array} \\right]\\left[ \\begin{array}{ccc}\\sqrt{10}&0\\\\0&0\\\\0&0\\end{array} \\right]\\left[ \\begin{array}{cc}\\frac{1}{\\sqrt{2}}&\\frac{1}{\\sqrt{2}}\\\\\\frac{1}{\\sqrt{2}}&-\\frac{1}{\\sqrt{2}}\\end{array} \\right]\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;可知：\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "\\Sigma = \\left[ \\begin{array}{ccc}\n",
    "\\sqrt{10} & 0 \\\\\n",
    "0 & 0 \\\\\n",
    "0 & 0\n",
    "\\end{array} \\right] \n",
    "= \\left[ \\begin{array}{c}\n",
    "1 \\\\ 0 \\\\ 0 \n",
    "\\end{array} \\right]\n",
    "[ \\sqrt{10}, \\ 0] =ab^T \\\\\n",
    "Ua = \\left[ \\begin{array}{c}\n",
    "\\frac{1}{\\sqrt{5}}\\\\\n",
    "\\frac{2}{\\sqrt{5}}\\\\ \n",
    "0 \n",
    "\\end{array} \\right], Vb = \\left[ \\begin{array}{c}\n",
    "\\sqrt{5} \\\\\n",
    "\\sqrt{5}\n",
    "\\end{array} \\right]\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;所以，$(Ua)(Vb)^T = \\left[ \n",
    "\\begin{array}{ccc}\n",
    "1 & 1 \\\\\n",
    "2 & 2 \\\\\n",
    "0 & 0\n",
    "\\end{array} \\right] =A$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 习题15.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;搜索中的点击数据记录用户搜索时提交的查询语句，点击的网页URL以及点击的次数构成一个二部图，其中一个结点集合$\\{q_i\\}$表示查询，另一个结点集合$\\{u_j\\}$表示URL，边表示点击关系，边上的权重表示点击次数。图15.2是一个简化的点击数据例。点击数据可以由矩阵表示，试对该矩阵进行奇异值分解，并解释得到的三个矩阵所表示的内容。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"images/15-2-Search-Click-Data.png\" style=\"zoom: 50%;\">\n",
    "<br/><div style=\"color:orange; border-bottom: 1px solid #d9d9d9;display: inline-block;color: #000;padding: 2px;\">图15.2 搜索点击数据例</div></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答思路：**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 根据二部图写出矩阵$A$\n",
    "2. 计算矩阵$A$的奇异值分解\n",
    "3. 解释奇异值分解得到的矩阵$U、\\Sigma、V$代表的内容"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答步骤：**   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第1步：根据二部图写出矩阵$A$**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A=\\left[ \\begin{array}{cccc}\n",
    "0 & 20 & 5 & 0 & 0 \\\\\n",
    "10 & 0 & 0 & 3 & 0 \\\\\n",
    "0 & 0 & 0 & 0 & 1 \\\\\n",
    "0 & 0 & 1 & 0 & 0 \n",
    "\\end{array} \\right]\n",
    "$$\n",
    "其中，行向量分别对应$(q_1, q_2, q_3, q_4)$，列向量分别对应$(u_1, u_2, u_3, u_4, u_5)$ 。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第2步：对矩阵 $A$ 进行奇异值分解;**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[ 1.   -0.    0.   -0.01]\n",
      " [ 0.    1.    0.   -0.  ]\n",
      " [ 0.    0.    1.    0.  ]\n",
      " [ 0.01  0.    0.    1.  ]]\n",
      "S= [20.62 10.44  1.    0.97]\n",
      "V= [[ 0.    0.96 -0.   -0.    0.29]\n",
      " [ 0.97 -0.   -0.   -0.24 -0.  ]\n",
      " [ 0.24  0.    0.    0.97  0.  ]\n",
      " [ 0.    0.29  0.    0.   -0.96]\n",
      " [ 0.    0.    1.    0.    0.  ]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.array([[0, 20, 5, 0, 0],\n",
    "              [10, 0, 0, 3, 0],\n",
    "              [0, 0, 0, 0, 1],\n",
    "              [0, 0, 1, 0, 0]]) \n",
    "\n",
    "# 调用numpy的svd方法\n",
    "U, S, V = np.linalg.svd(A)\n",
    "\n",
    "print(\"U=\", U)\n",
    "print(\"S=\", S)\n",
    "print(\"V=\", V.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**第3步：解释奇异值分解得到的矩阵$U$、$\\Sigma$、$V$代表的内容**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;根据题意，可知搜索中的点击数据记录用户搜索时提交的查询语句为$\\{q_i\\}$，点击的网页URL为$\\{u_j\\}$，点击的次数为矩阵的值。根据奇异值分解定理，可知：\n",
    "- $U$表示用户输入的查询语句与网页特征之间的关系矩阵\n",
    "- $V$表示网页与网页特征之间的相似性\n",
    "- $\\Sigma$表示用户输入的查询语句映射到网页的权重"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;举个例子：从矩阵$U$的第一列可以表示查询语句$q_1$与网页特征1的关系密切，通过矩阵$V$可以表示网页$u_2$与特征1的相似度最高，所以查询 $q_1$的时候点击网页$u_2$的次数是最多的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&emsp;&emsp;奇异值分解的另一个作用就是提取网页的特征来将用户输入与网页映射到一个低纬度空间中。通过第2步的计算，可以发现$\\sigma_1,\\sigma_2$的值比较大，计算截断的矩阵外积展开$A' = \\sigma_1u_1v_1^T+\\sigma_2u_2v_2^T$，可以发现其与矩阵$A$大体上是相等的。奇异值下降的速度越快，那么矩阵包含的更多信息就越集中分布在前面几个值比较大的特征上面。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考文献"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "【1】张丽静,刘白羽,申亚男. 实对称矩阵对角化教学的应用案例[C]. 大学数学, 2019, 35(2):116-121  \n",
    "【2】方阵的对角化：https://www.matongxue.com/parts/4628  "
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "c50c7698ff25eb1d7dce4c52e3d5cc14b3b23d6a97d4fb70aad8c815eba6f63e"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "273.188px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
